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A system of coupled kinetic transport equations for the Wigner distributions of a free variable 
mass Klein-Gordon field is derived. This set of equations is formally equivalent to the full wave 
equation for electromagnetic waves in nonlinear dispersive media, thus allowing for the description 
of broadband radiation-matter interactions and the associated instabilities. The standard results for 
the classical wave action are recovered in the short wavelength limit of the generalized Wigner-Moyal 
formalism for the wave equation. 
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The Wigner-Moyal formalism provides an alternative formulation for non-relativistic Quantum Mechanics, where 
^ ' the wave function is replaced by a quasi-particle distribution function, the Wigner distribution, in phase-space, which 
evolves in time according to a transport equation for the Wigner-Moyal distribution 1]. One of the main advantages 
of the Wigner-Moyal formalism is the possibility to describe quantum fields with non-trivial statistical properties |2i] . 

Since the evolution of electromagnetic waves can be described, in the paraxial approximation, by a Schrodinger 
like equation is not surprising that the Wigner-Moyal formalism has been used to describe waves propagating in 
an inhomogeneous, dispersive, anisotropic and slowly varying in time medium However, the propagation of 
^ ■ electromagnetic waves is, in general, a two mode problem as clearly stated by the second order time derivative in the 
' wave equation for electromagnetic waves in any medium Q . Previous works have only dealt with the single mode 
, problem, formally equivalent to propagating a nonlinear Schrodinger-like field, but this approach clearly breaks down 

■ when backscattering of the field is an important ingredient for the full dynamics of the system, as it usually occurs in 
laser/radiation-matter interactions at high intensities. 

In this paper, we build an alternative description for the propagation of electromagnetic waves in dispersive and 
\ nonlinear media, capable of capturing both forward and backscattering dynamics of the electromagnetic waves, by 
generalizing the work developed by J. Javanainen et. al |Q] to a variable mass problem. In the short wavelength 
approximation the results regarding the classical wave action (e. g. 01) are recovered. 

We first construct a 2 x 2 Wigner matrix on the basis of the Hamiltonian form of the Klein-Gordon equation of 
' a charged scalar particle field, introduced by Feshbach and Villars The diagonal elements describe forward and 
"j^ , backward photon densities and off-diagonal elements correspond to cross-densities in phase-space. In the corresponding 

■ quantum problem, the mass is assumed to be fixed |^; here a further generalization is performed in order to study the 
' variable mass problem, since the response of the medium (in our case, a cold plasma) exhibits spatial and temporal 

^ ■ dependencies. The motivation of this approach is to develop the relativistic phase-space description of the dynamics 
of Klein-Gordon particles analogously to the Hamiltonian description used earlier by Bialynichi-Birula, Gornicki and 
Rafelski for the Dirac particles 9]. The advantage of this method, compared with the manifestly covariant descriptions 
' of [To] , is that a single time parameter is used to describe the dynamics; as a result, the phase-space is a genuine 
dynamical system with its fully prescribed dynamics, and its conservation laws. 

In order for the discussion to be self-contained, a brief outline of the Feshbach and Villars Hamiltonian formulation 
for the Klein-Gordon equation is given in Sectional along with the definition of the 2x2 Wigner matrix. In Section 
mil we present the derivation of the generalized Boltzmann-Vlasov equation of motion for the Wigner matrix, and the 
set of coupled equations for the four real phase-space densities, which are related to the diagonal and cross-densities of 
the Klein-Gordon particles and their antiparticles, are derived. A physical interpretation of the phase-space densities 
is provided. In Section llVI the short wavelength approximation of the Boltzmann-Vlasov equation of motion for the 
Wigner matrix is performed, and the well known results of the classical wave action (e. g. (3, and references therein) 
are recovered. Finally, in Section we state our conclusions. 
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II. THE FESHBACK-VILLARS FORMALISM FOR ELECTROMAGNETIC WAVES AND THE 

WIGNER MATRIX 

The propagation of an electromagnetic wave in a dispersive medium is described by the Klein-Gordon equation 



e'^df<i>-e^c'^Vl<P + Lul{r,t)^ = (1) 

Here, $ is stands for one of the spatial components of the vector potential, in the Coulomb gauge, and c is the speed 
of light in vacuum. The CGS system of units is assumed throughout the paper, except where explicitly noted. We 
observe that ojp plays the role of the mass of the Klein-Gordon field. For a cold plasma ujp{r, t), is the electron plasma 

1 /2 

frequency, a function of time and space, given by u!p{r,t) — (^Ane'^ne/me'y) where e is the electron charge, me is 
the electron mass, ne($(r, i), r, i) is the electron density in the plasma, and 7($(r, t), r, t) is the relativistic Lorentz 
factor of the plasma electrons. The electron plasma frequency can be an arbitrary function of both time and space, 
and it may be a source of nonlinearity pl| . Equivalently, Eq.(^ is also the starting point for nonlinear optics |12| . 
which traditionally uses the index of refraction instead of LUp to describe the properties of the media. 

A distinguishing feature of this problem is the presence of a positive dimensionless parameter e. This parameter can 
be taken inversely proportional to the scale size of the spatial inhomogeneities, being small, but finite, for a slowing 
varying medium. Wc stress that our discussion is valid for £ = 1, which corresponds to the usual wave equation. The 
inclusion of this small parameter will allow us to easily perform the analysis in the short-wavelength limit in Section 

To find the Hamiltonian form of Eq. we follow the Feshbach-Villars @ description, defining: 
cj>,x=^('f±—dt<A (2) 



V2 V ^po 

where w^q is the background electron plasma frequency, independent of time and space, such that LUp(r,t) = 



~p(r,t). Eq. can then be written as: 



ijj. 



2 2 ~ 2 

iedt(l) = -^Vli^ + x) + :^{^ + x)+^po^ (3a) 

ZLOpO ■iljJpO 



iedtx 



2uj. 



'pO 



V2(0 + x) + 



2uj 



pO 



(3b) 



There are several ways of separating equation in two coupled equations, respecting definitions Equation Q 
denotes the correct expansion for the scalar potential, ujp, and guarantees, unlike the prescription proposed in 
that the short wavelength limit leads to the standard classical wave action conservation, as shown in Section IV. 
If we now introduce ^, the two-component vector potential, as 



^ = 

equation |(2Jl can be rewritten as 

(t3 -I- iT2) 



2uj. 



'pQ 



(4) 



(5) 



where 7i,n is the Hamiltonian matrix. In equation and henceforth, the notation: 
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T2 = 


i 










-1 



(6) 



is used to represent the Pauli matrices. 
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Equation (0) is a Schrodinger-like equation. Thus, the usual Wigner formahsm can be apphed to a 2 x 2 Wigner 
matrix, instead of using the Wigner formahsm on a scalar function. Defining the 2x2 Wigner matrix using the 
Feshbach-Villars representation 



VF(k,r,i)= (^^^ ^^dye^^vl;(r_,t)vl/(r+,t) 



(7) 



where 4* — 4*^3 denotes the Feshbach-Villars adjoint of the vector potential 4*, r± = r ± y/2, and M'^ the usual 
three dimensional Cartesian space . By construction, W is an hermitian matrix in the Feshbach-Villars sense, namely 
W = = W, and its explicit form can be expressed in terms of the field components and Xi yielding: 



W = 
where: 



1 

V 27r£ 



-w* 



dye'-r (r+,i)0(r_,i) 



(8) 



(9a) 



<t>x 



V27r£ 



dye'-0* (r+,t)x(r_,i) 



(9b) 



\2TTe 



dye'— X* (r+,t)x(r-,i) 



(9c) 



It is useful to calculate the equation governing the time evolution of 4'. Taking the conjugate transpose of equation 
Eq. jSj, and multiplying on the right by ra, we obtain: 



(t3 + iT2) 



2w, 



(10) 



where we used r/ — i G {0,1,2,3} and {Ti,Tj} = 2To5ij, i,j G {1,2,3}. Equation ((Tn|l demonstrates that the 
Hamiltonian matrix is not hermitian. As a consequence the transport equation for the Wigner matrix will have to 
include both commutators and anticommutators. 

III. GENERALIZED BOLTZMANN-VLASOV EQUATION OF MOTION 

The equation of motion for the Wigner matrix can be derived from the Schrodinger-type equations ((Sj and 
(|10|l for, respectively, 4* (r_, t) and 4* (r+, t), following a standard procedure 13] ■ 
From equations H5I7I10|I we obtain 



dtW{k, r, t) = 7^l + 7^2 - -ujpo [t3, W] 
where 



(11) 



dtW{k,r,t)^ (^^^ J^^<iye'^ [at*(r_,i)*(r+,t) + *(r_,t)ft*(r+,t)] 



(12a) 



n - — (— ' 

^ ujpQ \2'Ke 



dye' ^ S (t3 + iT2) 
[ 

*(r_,t) 



V2 

^*(r-,i) 



V2 



*(r+,i) 



4'(r+,t)- 
(t3 -I- ira) 



(12b) 



4 



2eujpQ \2Tre 



dye'— [*(r_,<)^'(r+,i)(r3 +ir2)c:;2^ 
- (rg +zT2)J)^_*(r_,i)W(r+,t) 



(12c) 



where oj^^ — ujp{r±,t). Note that [A, i?] represents the usual commutator between operator A and B, and that we 
also used the following identities: 



Vr±5(r±,i) = Vr<?(r±,t) = ±2Vy.g(r±,i) 



(13) 



where g{r, t) stands for a first order differentiable function of r, and an arbitrary function of t. 

The first term on the right hand side of Eq. (|11() . 7?,i,is the kinetic contribution for the equation of motion, and it 
is always present, independently of the potential il'p(r, t), whereas the second term, 1^2, accounts for the the potential 
contribution. After some lengthy calculations TZi, can also be written as: 



7^l = 

where 
Ho(k) = 



pQ 



{ra + tT2,W} ~ - [Ho(k)(T3 + ira), W] + Q 



(14) 



(15a) 



ujpo \2TTe J 



dye' 



(t3 +ir2)Vy«'(r_,t) • Vy*(r+,i)- 
Vy«'(r_, i) • VyW(r+, i)(r3 + ira) 



(15b) 



and {A, B} represents the usual anticommutator. Consequently, the equation of motion already denotes the effects 
of the non hermiticity of the Hamiltonian Tim, as previously noted. 

The form of equation (|15b|l suggests that a second order derivative in space, and a commutator must be involved; 
in fact, Q can be re-written as: 



Q = 



lec 



leading to 



{T3+iT2)^,W 



7^l 



c^k • V,. , ^ i 

{T3+tT2,W} 



2lj. 



pQ 



iHo(k)(T3+ZT2),W 



-7^o(k)(r3+zr2),iy 



(16) 



(17) 



with the generalized wave number k'^ = k" — e 



2Vi 
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2 = k2 

To calculate TZ2, we observe first that the inverse Fourier transform of the Wigner matrix ((7|) yields the vector 
potential matrix 5'(r+, i)\E'(r_, t), viz.: 

(18) 



*(r+,t)^'(r_,t) = / dk'e-*— M^(k',r,i) 
Combining Eq. fll^ and Eq. H12c|l leads to 



7^2 = 



2eujpo \2-Ke J 7r3 



dk' / dye'- 



~,2 



iy(k',r,t)(T3+ZT2)ci; 

(r3+ir2)VF(k',r,t)w2^ 



Since Coj = e^'ai ■'^"■cijp(r, t), and in the case of a sufficient regular function g [H 



I'y''' g{iy)F{k')dk' 



iyk' 
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(19) 



(20) 
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Thus, 7I2 can be expressed as: 

T^2 = / COS (- Vr • Vk [W^, (t3 + ZT2)] + sin ( - V r • V k { W^, (T3 + ir2)} (21) 



is: 



Combining Eqns. H17I21() in Eq. Hll|l . the generaUzed Boltzmann-Vlasov equation of motion for the Wigner matrix 
dtW{\^, r, t) + (i) - ^) i {r3 + ZT2, Ty(k, r, t)} 



+ - 
where: 



Ho(k) + C i [(t3 + zra), W{\^, v, t)] + ^ [T3, Ty(k, r, t)] = (22) 

Z £ 

k-Vr, 5= ^^^sin(-Vr- Vk), C= ^^^cos(-Vr- Vk) (23) 

In the argument of the trigonometric functions, the operator Vr acts on the electron plasma frequency, and the 
operator Vk acts on the components of the Wigner matrix Vl^(k, r,t). The matrix operator 7io(k) reduces, in the 
short wavelength approximation, to the free particle Hamiltonian Hq (k) . Both C and S can be expanded in a power 
series, containing only even powers of £ and, therefore, 5* and C in equations H23|) have regular limits when £ — > 0. 
The operators H23|l can also be cast in an integral form, which is more useful for numerical implementations: 



(7, 5/(k, r, i) = / dk'/C± (k - k', r, <)/(k', r, t) (24) 
where 

Equation H22|l describes the equation of motion for the Wigner matrix, and therefore it describes the transport of 
W in phase-space. Even though it is formally equivalent to the wave equation it is more useful to represent our 
field in terms of real Wigner distributions. Since the Pauli matrices form a basis for the 2x2 complex matrix space, 
W admits the expansion 

W = ^(tqWo + iriVKi - iT^W^ + T3T^3) (26) 

where all the Wi, i E {0,1,2,3} are real distributions. Since W is not Hermitian, an expansion in the form W = 
J2t=oWiTi, where aU the W.„ i e {0,1,2,3}, are real functions, is not possibe. Nevertheless, W is Hermitian in 
the Feshbach-Villars sense and the natural expansion with real coefficients is given by Eq. H26|l . It is easily seen 
that Wo = 1^00 - W^^, Wi = 2lm{W^^), W2 = 2Re{W^^) and 1^3 = W^^ + W^^. The physical meaning of 
W2(k, r, t) + WsCk, r, t) is of interest, because it is related to the Wigner function of the vector potential $. In fact, 
with 

VF$*(k,r,t)= (^^^ ^^dye^^$*(r+,<)ci>(r_,i) (27) 
where $* (r, t) stands for the complex conjugate of $ (r, t), then 

iy$$ - T4^2(k, r, t) + Wsik, r, t) = W^^ + W^^ + 2ReiW^^) (28) 

The vector potential is completely characterised by the distribution function W^^ , which can be calculated for different 
electromagnetic field configurations. Furthermore, there is a one-to-one correspondence between the distribution 
and the vector potential $, apart from a constant phase shift (pQ. Defining the inverse Fourier transform of as: 

•^ky(y,r,*)= / dke-'^[Ty2(k,r,i)+W^3(k,r,t)] (29) 



6 



$(r,t) is determined by: 



$(r,t) = 



(30) 



In Eq. (|3U|I (j)o is determined from the initial conditions, and the inverse Fourier transform of the Wigncr function is 
assumed to be different from at i = 0. 

The imaginary part of the cross-density W^^^i classical equivalent of the quantum 'Zitterbewegung' 

and it represents the interference between forward and backward photons. For Wigner matrices of a free forward 
plane wave or a free backward plane wave, of momentum kg, Wi is zero. 

As a simple illustration, let us consider a vector potential described by the superposition of two plane waves 
$(r, t) = $0 exp [i(ko • r — ujot)] + $i exp [i(ki • r — cuit)], propagating in an uniform medium with iip(r, t) = 0. For 
this field: 



2 \ i^po 



It) 



(31a) 



2 \ LUpo J 2 \ / 

The corresponding components of the Wigner matrix are: 

(^0 + ^po)^ 



W^^ik, r, t) 



2 -fl5{k - ko) 



{ujpO+ujo){ujpO+uji) ko + ki 
$o'I'i cos (77)0 k 



(31b) 



(32a) 



W,^{k, r, t) = -^K5[k ko) - - ki) + ^ 



1 - 



.a;o_wi 

WpO 



sin (ry) > 5 ( k 



ko 



cos (77) 
ki 



(32b) 



W'xx(k,r,i) 



(wo — ^po)' 
4^po 



$g<5(k-ko) 



(a;i - ujpof 
4wp0 



$?^(k- ki) 



(cjpo - UJq){UpO - 0^1 ) 



2'^po 



$o<i>i cos (?7)(5 k - 



ko + ki 



(32c) 



where 7y = Ak • r - A<^i, Ak = ki - ko, A^ 



^pO- 



— ciJo, <5(k) represents the Dirac delta distribution and ^0,1 — 
The presence of two photon beams associated with the two plane waves is easily seen. Furthermore 

there is a "third wave" which results from the interference between the two plane waves. In this case, Wi is different 
from zero, since it accounts for the interference between forward and backward photons. Wq is associated with the 
effective photon density, i.e., it is positive when the forward photon density is bigger than the backward photon 
density. 

Using the decomposition W — Re(iy) -|-iIm(T4^) in equation H22() . one can derive a set of two coupled real equations: 



dt2Re{W) + {b- S)]- {t3 + ZT2, 2Re(P^)} - - \no{k) + 6\\ [(rg + irj), n] Wi 



(33a) 



5t(TiVKi 



Ho(k) + C 



\ [(r3 + ZT2), 2Re(M^)] + [ra, 2Re(W)] = 
2 £ 



(33b) 
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where we have used {Ti,Tj} = 2TQSij, i,j G {1,2,3} and 2Im(Vr) = tiWi. From expansion H26() it is easily seen 
that 2Re(PF) = tqWq — iT2W2 + and, therefore, the four real phase-space densities obey the set of transport 

equations: 



dtWo + {D- S){W2 + W^a) = 



(34a) 



dtWi 



Wo(k) + C 



{W2 + W^) - -uJpqW2 



(34b) 



dtW2 -{D- S)Wo + - \no{k) + c]Wi + -ujpoWi ^ 
e L J £ 



(34c) 



dtWa + {D - S)Wa - - Wo(k) + C 



Wi=0 



(34d) 



The set of equations (|34|l is the main result of this paper. With the appropriate definitions for Wi, the set of 
equations H34|l is formally equivalent to equation QJ, and therefore opens the way to an appropriate description of 
the propagation of electromagnetic waves with arbitrary spectral content in nonlinear dispersive media [T3 | . 



IV. THE SHORT WAVELENGTH LIMIT 



We now examine the short wavelength limit of the Boltzmann-Vlasov equation of motion 122|l by examining the 
limit £ — > 0. In order to consider the short wavelength limit of H22|l we first observe that: 



VrW^(r,i) -> 



limC = 



lim k = k 



(35) 



Let the upperscript '■"^ represent any function in the short wavelength limit. In the limit £ — > 0, the set of equations 
(|33|l converges to a meaningful limit only if Wi faster than £ ^ 0; hence W'^^'^ must be real, i.e. W^'^'^ = 0. 
The same conclusions can also be obtained using the coarse-graining technique introduced by Shin and Rafelski . 
Using the results in H35|l and the previous considerations, Eqs. H33|) can be written, in the e — s- limit, as: 



-k- V, 



1 VrCi>l{r,t) 



~2( 



^{(r3+iT2),M/W} 



(36a) 



(t"3 + iT2) 



pO 



^pOT3 



,iy(o) 



= 



(36b) 



Our first goal now is to find a matrix operator that transforms the equation (|36b|l . into a necessary condition for 
W' — U^^W'-'^'^U to be diagonal. This will significantly reduce the amount of calculations needed to identify the form 
of the classical wave action, since the system will be reduced to two degrees of freedom. The matrix operator that 
performs this task is: 



(37) 



where C = h log 



a)(k,r,t) 



induced by U in Eq. (|36b|l . we obtain: 

Lu{k,r,t) [t3,W']^0 

where 



and with ^(k, r, t) — ^/k^c^ -I- lu^q + ujp{r, t). On performing the similarity transformation 



U-\Ts+iT2)U 



t^(k, r, t) 



(t3 + iT2) 



(38) 
(39) 
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From equation H38() we can see that W must be a linear combination of tq and T3 , since all the other Pauli matrices 
do not commute with T3. One can choose this linear combination to be parameterized by two real functions / and g, 
such that: 

W = l^ro + (40) 
Performing the transformation 137|l in Eq. H36a|l . with W given by Eq. (|4()|l . the diagonal part leads to 

aj(k,r,t) + .■ ^. • Vr/(k,r,t) - • Vk/(k,r,t) = (41a) 

c^k -> lVr(I'2(r,t) 
5tg(k, r, t) - ■ Vr5(k, r, t) + - " ' ' ■ Vkff(k, r, i) = (41b) 

Equations H41a() and (|41b|l clearly resemble the classical wave action conservation for forward and backward radia- 
tion, respectively |7i. However to confirm this fact, the relation between /, g and W^^ needs to be established. On 
one hand, we know that Re(iy) is l/2(roWo — iT2W2 + T3W3). On the other hand, in the short wavelength limit, 
Bje{W) is W'^°\ Thus, we conclude that: 

W^'^ = I (ro<) - zr.W-f + r^W^^) (42) 
Performing the transformation induced by U in Eq. H42|) we obtain: 

2W' = {f- g)To + (/ + g)T3 = U-\toW^°^ - ir^wlf^ + T^wi°^)U (43) 
which after some manipulation can be reduced to the set of algebraic equations: 

f-g^W^""^ (44a) 



/ + 5 = ^T^r + ^l^r (44b) 



- ^^i"' + (44c) 
where we have defined R = w(k, r, t)/tjpo- From Eqs. one infers: 



aj(k,r,t) (0) _ c2 



/ + .9 = ^^^^<^ = :^J(k,r,t) (45) 



where J(k, r,t) represents the classical wave action, as defined in ref. 0. Thus, in the short wavelength limit, 
Eqs. H41al41b|) state that the classical wave action, J(k, r, t), is conserved along the geometrical optics ray, indepen- 
dently of the adiabaticity condition [3| ■ 

To conclude the discussion of the short wavelength limit, we observe that when e in Eq. lO X is negligible, 
since reflected waves are exponentially small [l^ . Dropping the small component x Eq. ^ and eliminating the 
cjpo term by setting 0(r, t) = 0(r, t) exp {—iujpot/e), yields the paraxial wave equation |l8l |: 

2iLUpoedf4> = -e^c^ + Cj^^) (46) 

This equation is formally equivalent to the Schrodinger equation, and the Wigner-Moyal formalism can be followed 
in a straightforward way. This is the approach followed by Besieris and Tappert |3, 13 ■ expected, this path also 
leads to classical wave action conservation in the same limit we have explored here, but fails to capture the dynamics 
of the backscattered component of the radiation field. 



9 



V. CONCLUSIONS 

In order to establish a description of the propagation of electromagnetic waves in nonlinear dispersive media, 
including the two modes described by the wave equation, we have employed the Feshbach-Villars description of the 
Klein-Gordon field with a variable mass term, and the corresponding 2x2 Wigner matrix. The equations of motion for 
the Wigner matrix lead to a set of four coupled equations for four real phase-space densities, whose physical meanings 
are clearly identified. The general Boltzmann-Vlasov equations of motion account for the backscattered waves, and 
for the interference between forward and backward radiation. These equations replace the standard Wigner-Moyal 
description of classical fields Q, generalizing it to the two mode problem with variable mass terms, and are formally 
equivalent to the full wave equation. 

The short wavelength limit of the Boltzmann-Vlasov equations of motion was considered, and it was shown they 
reduce to uncoupled equations for the forward and backward radiation phase-space densities. In this limit, the 
standard results regarding wave action conservation, are recovered. 

The formalism introduced here provides a complete description for the electromagnetic field in nonlinear dispersive 
media, and opens the way to the description of broadband radiation driven parametric instabilities in nonlinear media. 
Applications of this formalism will be presented elsewhere. 
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